Learning Outcomes
In [1]:
import warnings
import numpy as np
import scipy as sp
import openpnm as op
%matplotlib inline
np.random.seed(10)
ws = op.Workspace()
ws.settings['loglevel'] = 40
np.set_printoptions(precision=4)
pn = op.network.Cubic(shape=[10, 10, 10], spacing=0.00006, name='net')
When performing transport simulations it is often useful to have 'boundary' pores attached to the surface(s) of the network where boundary conditions can be applied. When using the Cubic class, two methods are available for doing this: add_boundaries, which is specific for the Cubic class, and add_boundary_pores, which is a generic method that can also be used on other network types and which is inherited from GenericNetwork. The first method automatically adds boundaries to ALL six faces of the network and offsets them from the network by 1/2 of the value provided as the network spacing. The second method provides total control over which boundary pores are created and where they are positioned, but requires the user to specify to which pores the boundary pores should be attached to. Let's explore these two options:
In [2]:
pn.add_boundary_pores(labels=['top', 'bottom'])
Let's quickly visualize this network with the added boundaries:
In [3]:
#NBVAL_IGNORE_OUTPUT
fig = op.topotools.plot_connections(pn, c='r')
fig = op.topotools.plot_coordinates(pn, c='b', fig=fig)
fig.set_size_inches([10, 10])
OpenPNM uses a list-based data storage scheme for all properties, including topological connections. One of the benefits of this approach is that adding and removing pores and throats from the network is essentially as simple as adding or removing rows from the data arrays. The one exception to this 'simplicity' is that the 'throat.conns' array must be treated carefully when trimming pores, so OpenPNM provides the extend and trim functions for adding and removing, respectively. To demonstrate, let's reduce the coordination number of the network to create a more random structure:
In [4]:
Ts = np.random.rand(pn.Nt) < 0.1 # Create a mask with ~10% of throats labeled True
op.topotools.trim(network=pn, throats=Ts) # Use mask to indicate which throats to trim
When the trim function is called, it automatically checks the health of the network afterwards, so logger messages might appear on the command line if problems were found such as isolated clusters of pores or pores with no throats. This health check is performed by calling the Network's check_network_health method which returns a HealthDict containing the results of the checks:
In [5]:
a = pn.check_network_health()
print(a)
The HealthDict contains several lists including things like duplicate throats and isolated pores, but also a suggestion of which pores to trim to return the network to a healthy state. Also, the HealthDict has a health attribute that is False is any checks fail.
In [6]:
op.topotools.trim(network=pn, pores=a['trim_pores'])
Let's take another look at the network to see the trimmed pores and throats:
In [7]:
#NBVAL_IGNORE_OUTPUT
fig = op.topotools.plot_connections(pn, c='r')
fig = op.topotools.plot_coordinates(pn, c='b', fig=fig)
fig.set_size_inches([10, 10])
The boundary pores we've added to the network should be treated a little bit differently. Specifically, they should have no volume or length (as they are not physically representative of real pores). To do this, we create two separate Geometry objects, one for internal pores and one for the boundaries:
In [8]:
Ps = pn.pores('*boundary', mode='not')
Ts = pn.throats('*boundary', mode='not')
geom = op.geometry.StickAndBall(network=pn, pores=Ps, throats=Ts, name='intern')
Ps = pn.pores('*boundary')
Ts = pn.throats('*boundary')
boun = op.geometry.Boundary(network=pn, pores=Ps, throats=Ts, name='boun')
The StickAndBall class is preloaded with the pore-scale models to calculate all the necessary size information (pore diameter, pore.volume, throat lengths, throat.diameter, etc). The Boundary class is speciall and is only used for the boundary pores. In this class, geometrical properties are set to small fixed values such that they don't affect the simulation results.
In [9]:
air = op.phases.Air(network=pn)
water = op.phases.Water(network=pn)
water['throat.contact_angle'] = 110
water['throat.surface_tension'] = 0.072
In [10]:
from openpnm.phases import GenericPhase
class Oil(GenericPhase):
def __init__(self, **kwargs):
super().__init__(**kwargs)
self.add_model(propname='pore.viscosity',
model=op.models.misc.polynomial,
prop='pore.temperature',
a=[1.82082e-2, 6.51E-04, -3.48E-7, 1.11E-10])
self['pore.molecular_weight'] = 116 # g/mol
self.add_model commands within the __init__ section of the class definition. This means that when the class is instantiated, all the models are added to itself (i.e. self).**kwargs is a Python trick that captures all arguments in a dict called kwargs and passes them to another function that may need them. In this case they are passed to the __init__ method of Oil's parent by the super function. Specifically, things like name and network are expected.
In [11]:
oil = Oil(network=pn)
print(oil)
In the tutorial #2 we created two Physics object, one for each of the two Geometry objects used to handle the stratified layers. In this tutorial, the internal pores and the boundary pores each have their own Geometry, but there are two Phases, which also each require a unique Physics:
In [12]:
phys_water_internal = op.physics.GenericPhysics(network=pn, phase=water, geometry=geom)
phys_air_internal = op.physics.GenericPhysics(network=pn, phase=air, geometry=geom)
phys_water_boundary = op.physics.GenericPhysics(network=pn, phase=water, geometry=boun)
phys_air_boundary = op.physics.GenericPhysics(network=pn, phase=air, geometry=boun)
To reiterate, one Physics object is required for each Geometry AND each Phase, so the number can grow to become annoying very quickly Some useful tips for easing this situation are given below.
Perhaps the most distinguishing feature between pore-network modeling papers is the pore-scale physics models employed. Accordingly, OpenPNM was designed to allow for easy customization in this regard, so that you can create your own models to augment or replace the ones included in the OpenPNM models libraries. For demonstration, let's implement the capillary pressure model proposed by Mason and Morrow in 1994. They studied the entry pressure of non-wetting fluid into a throat formed by spheres, and found that the converging-diverging geometry increased the capillary pressure required to penetrate the throat. As a simple approximation they proposed $P_c = -2 \sigma \cdot cos(2/3 \theta) / R_t$
Pore-scale models are written as basic function definitions:
In [13]:
def mason_model(target, diameter='throat.diameter', theta='throat.contact_angle',
sigma='throat.surface_tension', f=0.6667):
proj = target.project
network = proj.network
phase = proj.find_phase(target)
Dt = network[diameter]
theta = phase[theta]
sigma = phase[sigma]
Pc = 4*sigma*np.cos(f*np.deg2rad(theta))/Dt
return Pc[phase.throats(target.name)]
Let's examine the components of above code:
target object as an argument. This indicates which object the results will be returned to. f value is a scale factor that is applied to the contact angle. Mason and Morrow suggested a value of 2/3 as a decent fit to the data, but we'll make this an adjustable parameter with 2/3 as the default.pore.diameter is actually a Geometry property, but it is retrieved via the network using the data exchange rules outlined in the second tutorial.target like a Physics object, that is a subset of the full domain. As such, the last line extracts values from the Pc array for the location of target and returns just the subset.sigma = 'throat.contact_angle'. This way the user can control where the contact angle could be stored on the target object.As mentioned above, the need to specify a separate Physics object for each Geometry and Phase can become tedious. It is possible to copy the pore-scale models assigned to one object onto another object. First, let's assign the models we need to phys_water_internal:
In [14]:
mod = op.models.physics.hydraulic_conductance.hagen_poiseuille
phys_water_internal.add_model(propname='throat.hydraulic_conductance',
model=mod)
In [15]:
phys_water_internal.add_model(propname='throat.entry_pressure',
model=mason_model)
Now make a copy of the models on phys_water_internal and apply it all the other water Physics objects:
In [16]:
phys_water_boundary.models = phys_water_internal.models
The only 'gotcha' with this approach is that each of the Physics objects must be regenerated in order to place numerical values for all the properties into the data arrays:
In [17]:
phys_water_boundary.regenerate_models()
phys_air_internal.regenerate_models()
phys_air_internal.regenerate_models()
The pore-scale models are stored in a ModelsDict object that is itself stored under the models attribute of each object. This arrangement is somewhat convoluted, but it enables integrated storage of models on the object's wo which they apply. The models on an object can be inspected with print(phys_water_internal), which shows a list of all the pore-scale properties that are computed by a model, and some information about the model's regeneration mode.
Each model in the ModelsDict can be individually inspected by accessing it using the dictionary key corresponding to pore-property that it calculates, i.e. print(phys_water_internal)['throat.capillary_pressure']). This shows a list of all the parameters associated with that model. It is possible to edit these parameters directly:
In [18]:
phys_water_internal.models['throat.entry_pressure']['f'] = 0.75 # Change value
phys_water_internal.regenerate_models() # Regenerate model with new 'f' value
In [19]:
inv = op.algorithms.Porosimetry(network=pn)
inv.setup(phase=water)
inv.set_inlets(pores=pn.pores(['top', 'bottom']))
inv.run()
'top' and 'bottom' using the pn.pores method. The algorithm applies to the entire network so the mapping of network pores to the algorithm pores is 1-to-1.run method automatically generates a list of 25 capillary pressure points to test, but you can also specify more pores, or which specific points to tests. See the methods documentation for the details.plot_drainage_curve. If you'd prefer a table of data for plotting in your software of choice you can use get_drainage_data which prints a table in the console.After running, the mip object possesses an array containing the pressure at which each pore and throat was invaded, stored as 'pore.inv_Pc' and 'throat.inv_Pc'. These arrays can be used to obtain a list of which pores and throats are invaded by water, using Boolean logic:
In [20]:
Pi = inv['pore.invasion_pressure'] < 5000
Ti = inv['throat.invasion_pressure'] < 5000
The resulting Boolean masks can be used to manually adjust the hydraulic conductivity of pores and throats based on their phase occupancy. The following lines set the water filled throats to near-zero conductivity for air flow:
In [21]:
Ts = phys_water_internal.map_throats(~Ti, origin=water)
phys_water_internal['throat.hydraulic_conductance'][Ts] = 1e-20
conduit_conductance. The term conduit refers to the path between two pores that includes 1/2 of each pores plus the connecting throat.
In [22]:
water_flow = op.algorithms.StokesFlow(network=pn, phase=water)
water_flow.set_value_BC(pores=pn.pores('left'), values=200000)
water_flow.set_value_BC(pores=pn.pores('right'), values=100000)
water_flow.run()
Q_partial, = water_flow.rate(pores=pn.pores('right'))
The relative permeability is the ratio of the water flow through the partially water saturated media versus through fully water saturated media; hence we need to find the absolute permeability of water. This can be accomplished by regenerating the phys_water_internal object, which will recalculate the 'throat.hydraulic_conductance' values and overwrite our manually entered near-zero values from the inv simulation using phys_water_internal.models.regenerate(). We can then re-use the water_flow algorithm:
In [23]:
phys_water_internal.regenerate_models()
water_flow.run()
Q_full, = water_flow.rate(pores=pn.pores('right'))
And finally, the relative permeability can be found from:
In [24]:
K_rel = Q_partial/Q_full
print(f"Relative permeability: {K_rel:.5f}")
Ti).Vp = geom['pore.volume'][Pi] and Vt = geom['throat.volume'][Ti].